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O ■ ABSTRACT 

' The Monge-Ampere-Kantorovich (MAK) reconstruction is tested against cosmological 

A-body simulations. Using only the present mass distribution sampled with particles, 
and the assumption of homogeneity of the primordial distribution, MAK recovers for 
each particle the non-linear displacement field between its present position and its 
Lagrangian position on a primordial uniform grid. 

To test the method, we examine a standard ACDM A-body simulation with Gaus- 
sian initial conditions and 6 models with non-Gaussian initial conditions: a \ 2 model, 
£N| . a model with primordial voids and four weakly non-Gaussian models. 

^ ' Our extensive analyses of the Gaussian simulation show that the level of accuracy 

of the reconstruction of the nonlinear displacement held achieved by MAK is unprece- 
dented, at scales as small as ~ 3ft. -1 Mpc. In particular, it captures in a nontrivial way 
the nonlinear contribution from gravitational instability, well beyond the Zel'dovich 
approximation. This is also confirmed by our analyses of the non-Gaussian samples. 
Applying the spherical collapse model to the probability distribution function of 
^j-^ | the divergence of the displacement field, we also show that from a well-reconstructed 

displacement field, such as that given by MAK, it is possible to accurately disentangle 
' dynamical contributions induced by gravitational clustering from possible initial non- 

Gaussianities, allowing one to efficiently test the non-Gaussian nature of the primordial 
Q ■ fluctuations. 

' In addition, we test successfully a simple application of MAK using the Zel'dovich 

c/3 , approximation to recover in real space the present-day peculiar velocity field on scales 

of 8 h~ l Mpc. 

Although non trivial observational issues yet remain to be addressed, our numer- 
ical investigations suggest that MAK reconstruction represents a very promising tool 
to be applied to three dimensional galaxy catalogs. 

Key words: cosmology: dark matter - large-scale structure of the Universe methods: 
A-body simulations - numerical 
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1 INTRODUCTION 

Recent data from the Cosmic Microwave Background 
(CMB), from the redshift surveys of the distribution of 
galaxies on large scales, and from weak lensing surveys on 
the projected cosmic density field have allowed not only a 
precise determination of the main cosmological parameters, 
but they have also confirmed the consistency and validity 
of the cold dark matter (CDM) paradigm for structure for- 
mation. The values of cosmologi cal parameters have been 
pinpointed by the WMAP data iSpergel et alJ 120031 from 
the power spectrum of the CMB temperature fluctuations, 
by a variety of results from the 2dFGRS and the SDSS such 
as the degree of anisotropy of the redshift space correlation 
of bright galaxies or their degree of clusteri ng (respectively, 
IPeacock et all 1200 it iTeemark et"aT]l2004aft . or by lensing 



surve ys (for example, Ivan Waerbeke et al] 120011 IPen et alJ 
12003ft . 

Some of these results depend on priors assumed for 
the primordial matter power spectrum (for example a "tilt" 
or a running index), and also on the nature of the dark 
matter, which affects the transfer function. However, recent 
measurements of both the present-day density correlations 
and the CMB fluctuations are consistent with the gravita- 
tional instability theory. The theory asserts that the struc- 
tures we observe today have evolved from infinitesimal, adi- 
abatic, Gaussian, scale-invariant primordial den sity fluctua- 
tions generated in minimal models of inflation <|Peiris et alJ 
l2003t ikomatsu et al J 120031: iTeemark et alJl2004rJ) . and ~\n- 
dicate that dark matter is "cold" at least at scales greater 
than about 100 kpc. 

Various details of the paradigm still need to be con- 
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strained. An important aspect is the shape of the proba- 
bility distribution function (PDF) of the primordial matter 
density field and in particular its possible deviations from 
Gaussianity on scales relevant to CMB and/or to large-scale 
structure (he reafter LSS) of the galaxy distribution (see 
|Peebleslll983l for early motivations). In fact, the simplest 
single-field inflationary models have been shown to be able 
to produce only a very small amount of non- Gaussianity 
belo w that detectabl e by current experimen ts (for exam- 
ple, ^feMacen^ IgOO^; [A^uwh^^^iyGQffil)- However, a 
variety of physically plausible models, including those where 
the seeds for structure form ation, at least partly , takes the 
form of topological defects iAvelino et alJ [l998 and refer- 
ences therein) or those associated to multifield inflation (e.g. 
Linde fc Mukhanovll997l:lPeebleJl999tlB"ernardeau fc Uzanl 



2002) can pro duce significant relic non-Gaussianity (see also 



Bartolo et aPl 20041) 



Recently, constraints on primordial CMB-scale non- 
Gaussianity have been obtained from hig her-order statis- 
tics o f temperature CMB maps (for example iKomatsu et alJ 
2003) assuming a parametrisation which bears the form of 
a quadratic nonlinearity in the primordial gravitational po- 
tential. Additional limits are provided by the abundance of 
massive clusters both today and at moderately high red shifts 
JChiu et alll998HKoyama et alJl999tlRobinson et alJ200ot 
iMatarrese et alJl2000|l 

Although current CMB maps and cluster abundances 
are consistent with primordial fluctuations being Gaussian, 
certain subtleties should be kept in mind. First, the sensitiv- 
ity of a given test in detecting primordial non-Gaussianity 
depends on the parametrisation chosen for the deviations 
from Gaussianity (e.g. whether one describes deviations 
in the gravitational pot ential or directly in the density: 
IVerde fc Heavens! l2001ah . Secondly, at a fixed comoving 
scale, Gaussianizing or non-Gaussianizing biases that can 
be induced by data reduction and also by gravitational 
clustering, if the scale is small, need to be carefully as- 
sessed and controlled. Finally, non-Gaussianity might be 
strongly scale-dependent as for models including topologi- 
cal defects or other kinds of phase transitions in the early 
Universe. Therefore, as is the case for the determination 
of cosmological parameters, the measurement of primordial 
non-Gaussianity will also benefit from combining different 
methods which rely on different assumptions, are sensitive 
to different shapes of non-Gaussianity and probe different 
scales. Assessing the possibilities of an original method in 
obtaining the statistics of the primordial density field espe- 
cially on comoving scales of ~ 3 to ;> 10 h~ x Mpc is the 
subject of this paper. 

A variety of methods which use LSS rather than CMB 
for the recovery of the statistics of the primordial den- 
sity field have been developed recently, some of which rely 
on measure ment of higher-order correlations of the galaxy 
distr ibution ([Ver de fc H eavensll2001ai iFeldman et al. 2001; 
IScoccimarro et alJ I2004T) and so me others on the sta tis- 
tics of the cosmic shear alone jTakada fc Jainl 120041) or 
in combination with the bright X-ray cluster abundances 
iAmara fc Refregierl I2004T) . The main difficulty with ex- 
tracting primordial non-Gaussianity out of the statistics of 
the present-day LSS resides in the contribution to non- 
Gaussianity at lower redshifts b y the nonlinear gravita - 
tional evolution of the density field llJuszkiewicz et alll993l) : 



gravity can generate skewness and higher-o rder moments 
in a random, initial ly Gaussian density field (Peebles 1980; 
iBouchet et"ali]ll992l) . 

An alternative approach for testing non-Gaussianity, 
consists of evolving the present time density field inferred 
from the galaxy catalogues back in time to higher redshifts, 
using reconstruction methods. Attempts to detect primor- 
dial non-Gaussianity applying reconstruction methods to 
galaxy surveys h ave been made, for example using the 1.2 
Jy IRAS s urvey (INusser et alJ | l995tl or IRAS PSC redshift 



catalogue iMonaco et alJ 120001 . Under the assumption of 



a linear bias the statistics of the primordia l density field 
were found to be consistent with Gaussianity iMonaco et alJ 
2000). However, most of these reconstruction methods re- 
quire a significant smoothing of the evolved density field 
to make it linear. This restriction calls for a Lagrangian 
based reconstruction algorithm which could be applicable 
at smaller scales. 

Variational methods have already been extensively used 
for the reconstruction of the primordial density and the 
present peculiar velocity fields, the determination of masses 
of galaxies and clusters and of the cosmological param- 
eters. In these methods one starts from the present po- 
sitions of the galaxies taken as mass tracers or of the 
dark matter particles (in the case of simulations) and finds 
their peculiar velocities and subsequently their full or- 
bits. These peculiar velocities (or orbits) ca n be obtained 
by m inimisin g the Euler-Lagrange action JPeebles! Il98gl : 
IShava et al]ll995l iNusser fc Branchinill2000l) . by using the 



perturbative Eulerian or 


1993; 


INusser et al. 


1995; 



ing an o ptimal assignment problem with stochastic al- 



iMonaco et al-llioOOt) or by solv- 
i proble m 

gorithms llCroft fc Gaztanagalll997l) or deterministic algo- 
rithms, such as the Mong e-Ampere-Kantorovitch (MAK) as- 
signment method jFrisch et alJl2002l: iMohavaee et alJl2003l : 
iBrenier et al.ll2003l . hereafter FMB1 which is used in this 
work. 

A difficulty with the reconstructions based on Pee- 
bles' action modeling is their lack of uniqueness: given a 
present distribution of mass, many different primordial den- 
sity fields and past histories, all physically plausible, can be 
reconstructed. All of these solutions correspond to station- 
ary points of the action. This can be a serious drawback 
when the primordial positions for several million particles 
have to be reconstructed, because the number of solutions 
can increase accordingly. An additional difficulty related to 
the lack of uniqueness of the solution is that variational- 
based methods cannot guarantee that the full possible solu- 
tion space has been explored. In addition, fast least-action- 
based algorithms which could be realistically applied to large 
datasets do not exist whereas fast MAK algorithms have 
been developed and are applied for the first time here to 
datasets of the order of 2 million particles. Furthermore, it 
has been shown that MAK reconstruction is a well-defined 
problem at large comoving scales (~ a few Mpc) where 
multi-steaming can be neglected and has a unique solution 
(see FMB and references therein). 

In this work, we test the MAK reconstruction method 
against iV-body simulations. We first focus on a simulation 
with Gaussian initial condition, namely standard ACDM 
model. We check the ability of MAK as a reconstructor 
of the non-linear displacement field, which is the difference 
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between the present position of a particle in the simula- 
tion and its initial position on a uniform grid. We extend 
the analyses to simulations with non-Gaussian initial con- 
ditions, and check whether MAK can detect the primor- 
dial non-Gaussianities. The non-Gaussian models examined 
in this work are a \ 2 model where the density field is the 
square of a Gaussian as a proxy for a strong non-Gaussianity 
which is scale-independent in our simulated volume, a pri- 
mordial voids model (PVM) to depict a scale-dependent 
non-Gaussianity which is significant at small scales but neg- 
ligible on large scales, and 4 models with quadratic correc- 
tions to the primordial gravitational potential (Q models) to 
represent a minute scale-independent non-Gaussianity con - 
sistent with the CMB (see for example lKomatsu et alJ 2003). 
We finally interpret some of the results of the analyses using 
the spherical collapse model. 

The paper is organised as follows. Section[5]summarises 
the features of the Gaussian and non-Gaussian models stud- 
ied here and describes the setup of the corresponding N- 
body simulations. An outline of our reconstruction method 
is provided in Section|H] Section|3]deals with the reconstruc- 
tion of the Gaussian model. We demonstrate that MAK 
is an extremely good reconstructor of the non-linear dis- 
placement field measured in the simulation, despite possi- 
ble limitations due to shell crossing. Section is devoted 
to the application of MAK to our series of plausible non- 
Gaussian models. Since MAK is a nonlinear reconstructor, 
we find that the primordial non-Gaussianities are mixed up 
with non-Gaussianities developed by gravitational instabil- 
ity. A simple but efficient way of disentangling these two 
non-Gaussianities is illustrated in Section |S| in terms of the 
spherical collapse model. In Section |7| we summarise the 
main results of this paper and discuss the potential extra 
complications that could arise in the application of MAK to 
real galaxy catalogues. 



2 SIMULATING MODELS OF THE 
PRIMORDIAL DENSITY FIELD 

We first describe the general features common to our 7 sim- 
ulations: a Gaussian ACDM reference model, a x 2 model, a 
primordial void model (PVM) and 4 models with different 
levels of quadratic non-Gaussianities in the primordial grav- 
itational potentials (Q models). We then describe the setup 
of the non- Gaussian initial conditions. 

2.1 Common simulation parameters 

We employ the public vers io ns of the iV-body codes 
HYDRA iCouchman et alJ Il995ft and GADGET 
iSpringel et al.l T^OCJlE to simulate collisionless struc- 
ture formation in the ACDM cosmo logy presently fa voured 
by the WMAP observations (ISpergel et alJ l2003f) . with 
parameters given in Tabled When constructing the initial 
conditions, we input primordial density field to the non- 
linear transformations when they are needed to generate 
the non-Gaussian field, and then smooth with the CDM 
transfer function of adiabatic fluctuations. 

The normalisation of the Gaussian case is set so that 
the value of the linearly extrapolated ag at z — is 0.9. 
This value provides only an approximate match to the 



z = abundance of clusters with optically derived masses 
jBahcall et alJ l2003f) which is then overesti mated, but it 
agree s with the WMAP preferred values llSpergel et alJ 
l2003t see also iTeemark et aljfeotMbl and references therein) . 
We tune the normalisation of the non-Gaussian models so 
that they produce the same abundance of massive clusters 
as the Gaussian simulation. In the \ 2 model (resp. PVM) 
we end up with a slight deficit (resp. excess) compared to 
the Gaussian case, but the quadratic series of models agrees 
very well with the abundance given by the Gaussian model. 

The starting redshift for all simulations is z = 70. The 
softening length is kept constant in comoving coordinates 
(set to one tenth of the mean interparticle separation). In 
what follows, we distinguish for clarity between the "initial" 
(or "primordial" ) and "unevolved" density fields which cor- 
respond to 2 = 70 and z — oo respectively. Particles are 
distributed on a regula r grid at z = oo, a nd then displaced 
with the usual scheme IZerdovichl lll970l) corresponding to 
the gravitational potential of the desired primordial density 
field realized on a 128 3 mesh (64 3 for the Q models). In the 
PVM model, additional displacement is imparted in regions 
covered by the primordial voids. 

In the next paragraph, we review the main features of 
the 3 classes of non-Gaussian models whose initial density 
fields we shall reconstruct in Section 

2.2 x 2 model 

The x 2 model with one degree of freedom is phenomenolog- 
ically attractive because it is strongly positively-skewed: it 
could be a valid description of a possible scale-dependent 
non-Gaussianity present on cluster scales. 

To realize the corresponding density field, 

Pinit(x) = ¥W(x), (I) 

we square on the 128 3 grid the scale-free n s = —2.4 Gaus- 
sian random field ipinit to obtain a n s ~ —1.6 x 2 density 
field pinit- This slope is slightly shallower than the theoret- 
ical value (—1.8) because of the incomplete mode coupling 
res ulting from squaring the field on the gr id (see for exam- 
ple |Whiti[l99l |Robinson We then apply 
the appropriate transfer function. Moments of the one-point 
PDF of a x 2 density field with power-law power spectrum 
(this holds on sufficiently small scales given the isocurva- 
ture transfer function) are theoretically independent of the 
smoothing length which is employed to compute th em, once 
the shape of the sm oothing window has been set fe.g. lPeeblesl 
Il999t IWhit3 n"9991. In practice when one uses a computa- 
tional mesh, this independence does not hold because of the 
incomplete mode coupling. 

Immediately after displacing particles from grid using 
the Zel'dovich approximation, we measure higher-order mo- 
ments of the resulting initial density field in the simulation 
using a cloud-in-cell (hereafter CIC) first order particle- 
to-grid assignment scheme, and smooth the resulting field 
with a top-hat kernel of radius 8/i _1 Mpc. As expected, the 
density is strongly positively skewed: we find D$ t g = 1.8 
and a kurtosis Dj .g = 9 . 5. (T he notations and defini- 
tions are those of jPeebles! [l999t Z> 3 ,8 = {{Sg/a 8 ) 3 } and 
P>4,s = {(5s/as) 4 } — 3 {(8g/ag) 2 } 2 .) These values have to 
be compared to the numerical values of 2.0 and 9.1 respec- 
tively found on a 128 3 , 200 ft" 1 Mpc grid realisation of x 2 
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smoothed with an 8 h" 1 Mpc radius. Here, we employ the 
X 2 case as a typical example of models with a strongly non- 
Gaussian initial density field but sti ll linear at the starting: 
redshift over all resolved scales (see iMoscardini et aDll99ll 
for a similar simulation). 



2.3 Primordial voids model 

The physical motivation for primordial voids resides in 
a possible first-order phase transition occurring in the 
early Universe in the framework of extended inflation 
<La fc Steinhardtl 11989ft . In some cases, incomplete nucle- 
ation of bubbles of true vacuum can result in the late-time 
persistence of a network of voids of astrophysically relevant 
radii which are completely empty of matter by the end of 
inflation/epoch of reheating. The radii of the voids can then 
increase in comoving coordinates onc e the Universe has en- 
tered the matter-dominated era fsee lBertschingeilll985l for 
the appropriate scalings). 

Specifically, we e mploy for PVM the parameters used by 
Griffi ths et al.l <l2003f) , and consider the voids as fully empty 
regions immediately surrounded by thin compensating shells 
of the matter which were swept up during their expansion. 
At z = 0, the cumulative number density of voids: 

N v {>r)=Ar a with a = -3 (2) 

is integrated from r m i n = 10 to r max = 25 h~ l Mpc. The 
index a = —3 is expected in viable models of the first-order 
phase transition th at we assume responsible for th e primor- 
dial voids (see, e.g. lOcchionero fc Amendolalll994) . The up- 
per cutoff r max approximate ly corresponds to the largest 
voids observed locally llEl-Ad fc Piranll200fi |PeebleJl20ol 
iHovle fc Vogelevll2004) . and r m i n is a possible lower cutoff 
depending on details of the phase transition. The normalis- 
ing constant A is then set such that to obtain a 40 percent 
volume fraction of the voids i n the present-day Un iverse. 

We refer the readers to iMathis et all J2004I) for the 
simple implementation of PVM in iV-body simulations. In 
short, we first displace the uniform grid of particles with the 
Zel'dovich prescription exactly as it is done for the standard 
Gaussian models. Then, before the simulation is started, we 
realize the "void network" . We compute a list of voids with 
centres randomly placed in the box and radii set according 
to equation (|5J (we ensure the voids are not initially overlap- 
ping). Next, we simply move any particle which is located 
inside one of these voids to its edge, along a straight line 
passing over the void centre and this particle's Lagrangian 
position. We assign the local shell velocity to the particle 
and then start the simulation. 

On 8 hT 1 Mpc scales, the skewness and kurtosis of 
the high-redshift density field measured immediately after 
constructing the initial conditions are Ds t $ = —1.68 and 
1)4,8 = 10.11. In fact, the void+Gaussian initial density 
field is both strongly negatively-skewed and non-linear on 
the scales of about 3 h~ Mpc. As such, it constitutes an 
important testbed for the possible reconstruction of the den- 
sity fields for which shell crossing has already occurred at 
z = 70. 



2.4 Quadratic contributions to primordial 
gravitational potentials 

The last four non-Gaussian cases which we study in this 
work, originate from the same generic model. The primor- 
dial gravitational potential includes a quadratic (also called 
"non- linear" ) contribution weighted by the parameter /nl 
in addition to the Gau ssian field $i n it (see for example 
IVerdefcHeave^l2001ah: 



$tot = $init + /nl (*Lt - (Sfnit)). (3) 



Equation © has become a convenient parametrisation 
of a small primordial non-Gaussianity t hat would be ex - 
pected in the simplest models of inflation (Maldaccna 2003). 
(A non-Gaussianity written as a quadratic correction to 
the primordial density field would be closer to a model 
with topological defects.) It has been found that higher- 
order correlation functions of the CMB temperature maps 
are significantly more efficient than those of LSS to con- 
stra in a non-Gaussian i ty such as that given b y equation 
© JVerde et alJl2000bUVerde fc HeavensfeoOlal) . The value 
of /nl is currently co nstrained to lie between -58 and 134 
ijKomatsu et, al.l|2003F) . Here, we simulate models with /nl 
of ±50 and ±100 (/nl = —100 is used to have a symmet- 
ric set). Positive (resp. negative) values of /nl result in a 
deficit (resp. an excess) of the abundance of massive clus- 
ters at 2 = compared to the Gaussian, /nl = 0, case. Note 
that definition implies that the level of non-Gaussianity 
in $tot depends on both /nl and ($ 2 nit ). Some caution is 
therefore necessary when applying the nonlinear transfor- 
mation to $init on the initial condition grid so that the 
non-Gaussianity effectively realized in the simulation on LSS 
scales corresponds to that probed in the CMB maps with the 
same /nl- This is done in the following four steps. First, we 
realize a Gaussian random potential field "l?mit with n s — — 3 
(corresponding to a density field with n s = —1). Second, we 
normalise it so that its extrapolation to the I ~ 10 scales 
probed by COBE jBunn et all 11995) with conversion of 
the resulting rms to CMB temperature fluctuations assum- 
ing adia batic modes corres ponds to the value measured by 
COBE l|Wright et alJll994h . Third, we realize equation Q 
using the normalised $i n it- Finally, we transform back $ to t 
to a density field and normalise it to the adequate value for 
the matter power spectrum normalization, as, at the start- 
ing redshift. Thereafter, the setup of the initial conditions 
is the same as that for the Gaussian simulation. 

Compared to the PVM and x 2 cases, the quadratic 
models considered here fall within CMB constraints and are 
significantly closer to Gaussian models. As a result, we ex- 
pect that distinguishing between these models and a Gaus- 
sian primordial density field by reconstruction techniques to 
be much more difficult for the Q models than for the x 2 or 
the PVM models. 

Next, we present the main aspects of the method em- 
ployed to reconstruct the initial displacement (and density) 
field from the z — outputs of our Gaussian and the 6 
non-Gaussian models. 
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3 RECONSTRUCTION OF DENSITY FIELDS 
AT HIGH REDSHIFTS 

In this section, we start by reviewing the main features of the 
MAK reconstruction method. Then, we detail the various 
particle samples which we have used for the reconstructions. 



3.1 Overview of the algorithm 

On sub-Zi^ 1 Mpc to h^ 1 Mpc scales, significant multi- 
streaming takes place in the low-redshift Universe, so that 
the pressurel ess Eulerian fluid approach to structure for- 
mation fails (|Peebles)ll987l) . Reconstruction is a well-posed 
problem insofar as significant multistreaming has not oc- 
curred. Our main hypotheses are that the mapping q — > x 
between the Lagrangian coordinate, q at iinit, and the Eu- 
lerian coordinate of a particle, x(q, tanal)i can be written 
as the gradient of a convex potential "l?(q, i). The convexity 
guarantees that the mapping is one-to-one, and therefore the 
existence of the reciprocal map x — > q. In addition, these 
hypotheses imply that reciprocal mapping also derives from 
a convex potential 0(x, i), and t hat O and $ a re related by 
the Legendre-Fenchel transform iArnoldlll978r . The goal of 
the algorithm is to obtain the mapping O. This is addressed 
in the first paragraph that follows. The details of the com- 
putation of the reconstructed peculiar velocities is discussed 
in the second paragraph that follows. 



3.1.1 From cosmological reconstruction to the assignment 
problem 

Substituting the inverse Lagrangian map, q = V x @(x), in 
the equation of conservation of mass, p(x)dx = p dq yields 
the so-called Monge-Ampere (MA) equation: 



det 



<9 2 9(x) 



dxidxi 



p 



(4) 



where p is the unevolved (uniform) Lagrangian density and 
p(x) the evolved Eulerian density. The linearisation of equa- 
tion @ yields the Poisson equation. Recently, it has been 
shown that the map that is solution to the MA equation is 
the unique solution to the Monge-Kantorovich optimisation 
problem (hence "MAK"), in which one seeks the map x — > q 
which minimises the quadratic cost function: 



1= / p|x-q| 2 d 3 g = / p(x)|x-q| 2 d 3 x 



(•>) 



fsee lBenamou fc Brenierb oOO. FMB, and references therein 
for the detailed proofs). In practice, when the method is ap- 
plied to discrete points such as galaxies or iV-body particles 
at z = 0, one discretizes the cost function © into: 



mm 

m 



) -Xil 



(6) 



which is known as the assignment problem. We stress here 
for clarity that both the set of initial and final positions at 
Zinit and Zfl na i are known, but the correspondence between 
the two is not. The set of final positions is measured from the 
observations or the simulations, that of the initial positions 
is produced using the fact that the unevolved distribution 
is uniform: in our case, it corresponds to the nodes of a 



regular mesh. Given iV initial and N final entries, the aim is 
to find the permutation which minimises the quadratic cost 
function 1 . 

The simplest algorithm which would solve the assign- 
ment problem (|SJ would clearly have a factorial complexity: 
one needs to search among N\ possible permutations for the 
one which has the minimum cost. However, advanced as- 
signment algorithms exist which reduce the complexity of 
the problem from factorial to polynomial. The latest algo- 
rithm developed by M. Henon and used in this work, which 
is a cosmological adaptation of Bertsekas' auction algorithm, 
scales approximately as A" 2,5 (for relevant details see for ex- 



ample lHenonl ll992 , 



Bertsekas 1998). 



3.1.2 Reconstruction of the peculiar velocity field 

Once we have obtained the optimal one-to-one mapping be- 
tween the present-day and its Lagrangian position for every 
particle, we can evaluate the present, reconstructed peculiar 
velocities for these particles using the Zel'dovich approxima- 
tion. From there, their reconstructed positions can be inter- 
polated at any desired redshift between Zjnit and Zfinal- This 
works as follows. In the Zel'dovich approximation, peculiar 
velocities are given by : 



x = /(fi) H(t) (x - q) 



(7) 



where f(fl) = dlnD/dlna is the dimensionless linear growth 
rate, D(t) is the amplitude of the growing mode today, a is 
the cosmic scale factor and H(t) is the value of the Hubble 
parameter. We note Do and xo the growth factor and the 
known present-day position of a particle respectively. If we 
now assume the validity of the Zel'dovich approximation 
throughout, we can integrate equation Q backwards from 
today (recall that q is known from the MAK mapping) and 
obtain the position at any z using: 



*(*) = q + -FT (xo - q) 
Do 



(8) 



Equations Q and JHJ can be regarded as a "backwards" 
Zel'dovich scheme going from z = to high redshifts to dif- 
ferentiate it clearly from the "forward" Zel'dovich scheme 
which is employed to go from z — oo to lower redshifts (e.g. 
z = 70 in our case in order to set the initial conditions 
of the simulations). As we look for the map satisfying the 
MA equation in the following sections, we start from the 
Zfinal = distribution of particles in the simulations and re- 
late it to an unevolved position which is a node of a regular 
lattice. Before we turn to the detailed analysis of the recon- 
structed high-redshift fields, we define in the next paragraph 
the series of particle samples that we have used for our re- 
constructions. 



3.2 Defining particle samples for reconstruction 

For each simulation, we have traced back a set Si of 64 3 par- 
ticles. These particles are regularly selected at z = oo from 
the uniform 128 3 grid of particles, then located at z = 
in the full, 128 3 simulation output, and finally used for the 



1 Because the simulation boxes are periodic, optimal assignment 
is found after taking his periodicity into account. 
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reconstruction. In the case of the Q series of models with 
64 3 particles, there is no resampling needed. The mean in- 
terparticle distance in Si is 3.1 h~ x Mpc. One motivation 
for this was to reduce the computational cost of our 7 re- 
constructions. 

In the Gaussian and the \ 2 cases, we have also recon- 
structed the full set So of 128 3 particles, to help our theoret- 
ical interpretation of MAK and to verify that our results are 
not affected by the under-sampling performed in Si . Finally, 
for the Gaussian, the x 2 and the PVM models, we have re- 
constructed a smaller, 64 3 particle "sub-box" extracted from 
the 200 /i -1 Mpc box, we call it "dense samples" S2; it is 
obtained as follows. First, we have extracted a 100 ft -1 Mpc 
-size square box of 64 3 particles from the original Lagrangian 
grid, and located them in the final simulation output. Then, 
we have reconstructed their final positions back in time. We 
will use this sample together with the sample So as we dis- 
cuss issues of resolution. The summary of the number of 
particles, boxsize and mean interparticle separation I of our 
various particle samples is given in Tabled 

The reconstruction of the 128 3 particles of the Gaussian 
and x 2 sample So is unprecedented and breaks a computa- 
tional barrier for cosmological reconstruction schemes. 

In Fig. we show for this simulation the scatter plot 
of the MAK-reconstructed Lagrangian positions of particles 
against their true Lagrangian positions, as well as the his- 
togram of the differences between true and reconstructed po- 
sitions. This very simple preliminary analysis already shows 
the success of MAK reconstruction, at least on large scales. 
However, we see that "exact" reconstruction is achieved for 
only a small fraction of particles. The collapse of highly non- 
linear structures indeed brings about shell crossings and lo- 
cal mixing. MAK is therefore unable to distinguish between 
particles belonging to the same condensed object as will be 
illustrated in detail later. In fact, we shall see that this does 
not pose an obstacle to the recovery of the primordial den- 
sity field. 



4 EFFICIENCY OF MAK AS A NONLINEAR 
RECONSTRUCTOR 

This section focuses on the Gaussian model. Firstly, 
we present detailed comparisons between the MAK- 
reconstructed displacement field with the non-linear dis- 
placement field (here after SIM NL) which it is designed to 
recover. Secondly, we compare the reconstructed and sim- 
ulated nonlinear "density fields" that will be given by mi- 
nus the divergence of the displacement field. We denote the 
latter by v and for this reason on occasions we refer to it 
as the "velocity field". In both cases, comparisons with the 
Zel'dovich, initial conditions fields (hereafter SIM IC) are 
also given. Thirdly, we compare the present velocities found 
by MAK to the true peculiar velocities given by the simula- 
tion: the agreement on large scales is a notable achievement 
of our reconstruction method. In the last paragraph, we dis- 
cuss potential resolution issues. 

Except for Section ^, 3l which concerns the present pecu- 
liar velocity field, the analyses are performed on the displace- 
ment field which is a Lagrangian quantity, always estimated 
as a function of Lagrangian coordinate, q, on the primor- 
dial grid of initial unperturbed positions of the particles. In 



Table 1. Parameters of the various particle samples (full box, 
sparse samples and dense sample, respectively So, Si and S2) of 
our 7 models used to test reconstruction. Note that although we 
have called "full box" the particle samples of the Q models which 
only have 64 3 particles in total, their mean interparticle separa- 
tion, I, makes them more similar to the "sparse sample" cases. As 
a result we also call them Si . Also given is the fraction of parti- 
cles that the MAK reconstruction scheme reassigns to their true 
Lagrangian position on the uniform grid. We call these particles 
"ideally reconstructed" . 



Sample N 


Box 

[ h- 1 Mpc ] 


I 

[ h- 1 Mpc ] 


ideal 


f2m = 0.3 , h = 0.65 , (Tg = 


Gaussian 
= 0.99 , n s = 1 , M = 3.2 x 10 11 


h- l M e 


full box (So) 128 3 
sparse sample (Si) 64 3 
dense sample (S2) 64 3 


200 
200 
100 


1.5 

3 

1.5 


18% 
37% 
17% 


Q m = 0.2, h = 0.7, cr 8 = 


x 2 

0.61 , n s = -2.4 , M = 2.1 x 1C 


i 11 /i- 1 M0 


full box (So) 128 3 
sparse sample (Si) 64 3 
dense sample (S2) 64 3 


200 
200 
100 


1.5 

3 

1.5 


31% 
60% 
33% 


Q m = 0.3, h = 0.7, a 8 = 


PVM 
0.9, n s = 1, M 


= 3.2 x Wall- 


■'Mg 


sparse sample (Si) 64 3 
dense sample (S2) 64 3 


200 
100 


's 

1.5 


28% 
10% 


Q m = 0.2, h = 0.7, cr 8 = 


Q models 
0.9, rig = 1, Af 


= 2.5 x W 1L h- 


■'M 


Q-100 


full box (Si) 64 3 


200 


3 


43% 


Q-50 


full box (Si) 64 3 


200 


3 


43% 


Q+50 


full box (Si) 64 3 


200 


3 


42% 


Q+100 


full box (Si) 64 3 


200 


3 


41% 



all cases, the measurements use the actual displacement be- 
tween z — 00 (the unperturbed grid) and z = 0. For SIM IC, 
this requires a scaling of the linear displacement field using 
the Zel'dovich approximation. 

Note furthermore that all the analyses in this section 
and in forthcoming ones are performed with additional top 
hat smoothing applied to the fields with spherical windows 
of various radii (3 ft -1 , 8 h~ x Mpc, and sometimes 12 h~ x 
Mpc). It is important to emphasize that top hat smooth- 
ing with a spherical window of radius R is approximately 
equivalent to Gaussian smoothing with a window of radius 
R/y/E. This should be kept in mind when comparing with 
previous works which usually use Gaussian smoothing (e.g. 
Monaco et al. 2000). 
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Figure 1. Performances of MAK reconstruction, left panel: Scatter plot of the MAK-reconstructed initial coordinates of a particle 
versus its initial, post-Zel'dovich coordinate for the 128 3 particles (So sample) in the Gaussian simulation. We use a "quasi-periodic 
projection" (QP) coordinate q = (q\ + q2\f 7 2 + g3v3)/(l + v2 + \/3), with q\ , etc... 6 [0,1] where 1 corresponds to the box size. These 
QP coordinates guarantee a 1-to-l correspondence between the q-values and the points on a regular grid. To clearly show where the 
majority of particles reside on the diagram, we show the decimal logarithm of 1 plus the local density of points (the resolution is 1/256 in 
QP coordinate, and the Lagrangian mesh spacing is 1/128). Note that most of the particles lie almost perfectly on the diagonal already 
illustrating the success of MAK. right panel: Histogram of the distribution of the distances, expressed in units of mesh size, between the 
true Lagrangian, discrete, positions of each particle on the uniform grid and the position assigned to it by the MAK reconstruction. The 
first bin corresponds to what we call "ideal" reconstruction, as shown in Table 1. If one adds up the percentages obtained in the first 
and the second bin, one should obtain a good estimate of the percentage of "ideally" reconstructed positions for the sparse sample Si, 
as is the case (see Table 1). 



4.1 Comparison of the displacement fields 

In the present and the next paragraph we analyse the Si 
sample of the Gaussian model to show that, as expected 
from the description of the algorithm in Section [3] MAK 
is an excellent tracer of the non-linear displacement field. 
We recall that the SIM NL displacement is simply mea- 
sured in the simulation by joining a particle's position on 
the Lagrangian grid to that particle's position in the z — 
simulation output. 

Fig. H compares, for SIM IC, MAK and SIM NL, the 
z-components of the displacement field projected along the z 
axis from a 20 Mpc thick slice. The measurements have 
been performed using a nearest grid point (NGP) interpola- 
tion of the particles to a 64 3 grid, followed by smoothing with 
a top-hat of radius 3 Mpc . There is a very good agree- 
ment between the three panels of the figure. To have addi- 
tional insight, Figure [!|] displays v z measured along a line. 
The two panels correspond to smoothing scales 3 h~ x Mpc 
and 8 Mpc . This figure highlights the subtle differences 
between the three panels of Fig. [5] In particular we see that 
the agreement between SIM NL and MAK is better than 
between SIM NL and SIM IC or between MAK and SIM 
IC. Note that even at mildly nonlinear scales, 3 hT 1 Mpc, 



the agreement between MAK and SIM NL is surprisingly 
good. 

These results are confirmed by Fig. [I] which examines 
the PDFs of the n z component of the displacement field, for 
SIM NL, MAK and SIM IC, at 3 ft" 1 Mpc : the curves cor- 
responding to MAK and SIM NL are nearly superimposed. 
They present a slight skewness which is due in the SIM NL 
case to non-Gaussianity brought by gravitational clustering. 
We already see that the latter feature is entirely captured 
by MAK. 

So far we have compared single components of the dis- 
placement field: we now address the differences and/or simi- 
larities in amplitude and in direction between the three dis- 
placement fields as a function of the underlying over-density. 

Figure displays the measured distribution of the ra- 
tios vsim ic/«mak , !)simnl/"mak and vsim ic / vbim nl as a 
function of initial density estimated as minus the divergence 
of the initial displacement field. The upper and lower rows of 
panels correspond to a smoothing scale of 3 and 8 ft -1 Mpc 
respectively, prior to the calculation of their ratios. Again, 
the amplitudes of the MAK and SIM NL displacements 
agree well. The match is slightly less good for initially over- 
dense regions, as expected. Indeed, these regions might have 
experienced shell-crossing during dynamical evolution. Since 
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Figure 2. Projected v z displacement field taken from a 20 h~ 1 Mpc -thick slice normal to the z-direction, for the Gaussian model. 
The resolution before projection is ~ 3 hT 1 Mpc in each direction. The left panel is the Zel'dovich initial conditions (SIM IC), the 
central panel is the MAK reconstruction and the right panel is the nonlinear displacement SIM NL. (Note that we will keep this tryptich 
structure and the corresponding SIM IC, MAK and SIM NL association for all remaining colour images.) The colour scale shown at the 
top is common to all plots and expressed in multiples of the rms of v z of each field. The largest and broadest features repeat well in all 
three maps, but the MAK and SIM NL maps also match in detail. 



SIM NL and MAK agree rather well, we notice similar dis- 
crepancy between SIM NL versus SIM IC and MAK versus 
SIM IC. 



Fig. is the same as Fig. but for the angle between 
the smoothed components of the displacement fields. The 
measured angle is always positive, but we made for clarity a 
reflection at the horizontal axis. To help understanding this 
figure one can assume for instance the SIM IC displacement 
direction to be chosen at random within a cone of angular 
size 10 (respectively 20) degrees centred on the MAK/SIM 
NL displacement direction. In this case, the observed distri- 
bution of angles would be bimodal and would peak at 6.7 
(respectively 13.3) degrees. 



In regions of low to average densities, the quality of 
the MAK reconstruction of the directions of the SIM NL 
displacement is striking. Above the mean density, the scat- 
ter starts to increase as expected again because these zones 
roughly correspond to virialized regions at z = 0, where 
reconstruction is more difficult, but the bulk of the PDF re- 
mains within a ± 2 degrees offset. Note also that in the most 
strongly over-dense regions the displacement amplitudes are 
typically smaller than in regions of lower density, so that the 
consequences of significant errors in direction are not dra- 
matic, simply because the magnitude of the displacement is 
small. On the other hand, comparing the MAK/SIM NL to 
SIM IC displacements, we find that the agreement is much 
less good, the most likely difference between the two direc- 
tions being larger than ~ 2.5 degrees at the smoothing scales 
considered here. 



4.2 Comparison of the "density fields" 

We estimate the density contrast, S, using the linear-regime 
relation <5(x) oc — V 9 ■ v. For the SIM IC sample this is a fair 
estimate of the density field at this redshift. For the SIM NL 
and MAK samples, this basically gives an estimate of the 
primordial density by simple application of the Zel'dovich 
approximation. 

In practice, to estimate V 9 • v, we first assign the dis- 
placement field of the particle set Si to a 64 3 grid using 
nearest grid point interpolation (hereafter NGP) and then 
Fourier transform to take the gradient. As a result of tak- 
ing the derivative, the V 9 • v field is intrinsically more noisy 
than the displacement field. Fig. [7] is the pendant of Fig. [5] 
but for V 9 • v. It is supplemented with Fig. |S] which shows 
the 3- and 8 h~ l Mpc smoothed divergence field computed 
along a line parallel to the x direction cutting through the 
centre of the box. 

The MAK density field clearly reproduces the SIM NL 
field very well at all values of V 9 ■ v. On the contrary, the 
MAK field only reproduces well the SIM IC field in regions of 
low to average densities. High-density (y ~ 1.5 — 2) patches 
are much more extended in MAK and in SIM NL than in 
SIM IC, and the v ~ 2.5 to 3 levels of the SIM IC field are 
barely reached by MAK or SIM NL. 

This can be interpreted as follows. The fact that the 
simulated and the reconstructed displacement fields agree 
very well in the under-dense regions in all the cases is not 
surprising: indeed, the Zel'dovich mapping, hence SIM IC, 
is known, at l east to some approximatio ns, to work well in 
these regions iSahni fc S handarin 1996) and MAK recon- 
structs perfectly the Zel'dovich dynamics prior to the shell 
crossing. Therefore, in the under-dense regions, agreement 
among SIM IC, SIM NL and MAK is expected at least qual- 
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Figure 3. Profiles of the initial z-component of the displacement fields v z in the Gaussian model, measured along a line parallel to the 
x-direction and located in the middle of the box. The resolution is 3 (resp. 8) h" 1 Mpc in the upper (resp. lower panel). The Zel'dovich 
SIM IC, MAK reconstruction and SIM NL profiles are shown with dotted, dashed and dash-dotted lines respectively. The horizontal 
dotted lines give ±cr of the SIM IC variations (the rms is measured from the profiles). Note how well MAK reproduces the SIM NL 
displacement field, even at 3 h" 1 Mpc resolution. 
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relative likelihood given density [ arbitrary units ] 




0.000 0,196 0.390 D.5B6 0.760 0,976 1.170 




-a -z -l o i e a "3—3—10 i e a —a —3 — i o i e a 

- div T/tFeorv [ unit Ems ] in SIM IC - dlv v/fro*, [ unit rma ] in Sill IC - div t/ct^t [ unit rma ] in SIM IC 

Gaussian 



Figure 5. Left panels show the conditional PDF of the local ratio between the amplitudes of the Zel'dovich SIM IC and of the MAK- 
reconstructed displacements, as a function of the 3 h~ 1 Mpc local over-density in the initial conditions SIM IC. The plot corresponds 
to the Gaussian model. The top (resp. bottom) row is for a 3 (resp. 8) h -1 Mpc smoothing length. Middle panels: same as the left 
panels, but for the ratio between the amplitudes of the non-linear displacement field SIM NL and the reconstructed displacement field. 
Right panels: same as the left panels, but showing the ratio between the amplitudes of the Zel'dovich, SIM IC displacement field and 
the non-linear displacement field SIM NL. 



itatively. In detail, however, extra complications arise due to 
nonlinear dynamics which will be discussed below in more 
detail (see also Fig. II IB . 

In the over-dense regions, the interpretation of the re- 
sults is more complex: the large patches with roughly con- 
stant density in SIM NL and MAK correspond to collapsed 
objects at z = 0, as explicitly illustrated by Fig. [5] where 
we analyse a typical cluster at full resolution, i.e. extracted 
from the So sample. Using Fig. as a guideline, let us con- 
sider a Lagrangian region that will collapse to an object of 
negligible size, as illustrated by the upper left panel. It is 
trivial to realise that the displacement from the initial to 
the final stages is purely radial as shown in the left and 
right bottom panels for MAK and SIM NL respectively, and 
therefore that minus its divergence is equal to 3, which is 
the upper bound expected for MAK reconstruction and for 
the simulation (more approximately in the latter case). 

We thus see that the left and right bottom panels 
of Fig. EH are extremely similar. They mostly differ in the 
boundary of the Lagrangian region occupied by the cluster: 
this means that the main difference between MAK and SIM 



NL relies on particles which have been assigned to the wrong 
collapsed objects. Particles that are not in the intersection 
of the two patches of the left and right bottom panels are 
clearly assigned wrongly by MAK and will give totally wrong 
displacement fields. We, however, notice that the boundaries 
of these "peak patches" are nearly the same, which explains 
why MAK performs so well. At this point, it is important 
to notice that exact recovery of initial positions of particles 
belonging to the same collapsed object is not crucial and 
the quality of the reconstruction of the displacement field is 
only marginally affected if the positions of the particles end- 
ing up in the cluster are swapped by the MAK assignment. 
This demonstrates that there is no need to worry about the 
occasional low level of "ideal" reconstruction presented in 
Table 1. Of course, the linear displacement field is not ra- 
dial as shown in the upper right panel of Fig. |U] because 
it accounts for the details in the distribution of the initial 
fluctuations that are in the Lagrangian volume: in a hier- 
archical model such as CDM, smaller structures form first 
which then merge and create larger structures. In our pic- 
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Figure 6. Same as Fig.|^]but for the angle [in degrees] between the displacement fields indicated by the vertical label of each plot, as a 
function of the over-density smoothed on a 3 h~ 1 Mpc scale. 



ture, this would mean successive episodes of approximately 
locally radial displacements. 

Clearly, Fig. |7| already emphasises the strengths and 
weaknesses of MAK: it seems to reconstruct the nonlinear 
displacement field and its divergence very well. However, the 
real goal is the recovery of the primordial density field: in 
this case, MAK performs qualitatively well in the under- 
dense regions but fails to do as well in the high-density 
regions due to the presence of collapsed structures, as ex- 
plained above. 

To illustrate these results in a more quantitative way, 
Fig. 1 1 Ul examines the PDF of the divergence of the displace- 
ment. The upper and lower panels of the figure correspond 
to the field smoothed here within a radius of ~ 3 ft -1 Mpc 
and ~ 8 ft -1 Mpc respectively. Again, the PDF of the MAK 
density field reproduces very well that of the SIM NL density 
field, including its positive skewness and its abrupt cutoff at 
large v which corresponds to the present value of — V 9 -v = 3 
as discussed above. 

Fig. Illl goes into deeper detail by examining the bivari- 
ate distributions: it represents scatter plots of the divergence 
of the displacement field for MAK versus SIM IC (left pan- 
els), MAK versus SIM NL (middle panels) and SIM NL ver- 
sus SIM IC (right panels) computed for the fields smoothed 
at 3 and 12 h~ Mpc (upper and lower rows). The banana- 



shape observed in the right panels is due to nonlinear cluster- 
ing and can be predicted from perturbation theory, similarly 
to what has p reviously been done for the divergence of th e 
velocity field feernardeaulll994bl iBernardeau et al"lll999T) . 
A calculation of — V 9 • v relying on the spherical collapse 
model is given in Section |5] and displayed on the bottom 
right panel. If one neglects the cutoff expected for SIM NL at 
— V 9 • v ~ 3 the nonlinear gravitational effects increase the 
value of — V 9 ■ v in the under-dense regions and increase it 
in the over-dense regions compared to the linear predictions. 
We notice here that the qualitative results stated previously 
about under-dense regions need to be quantified: even in 
the low density regime the recovery of the initial density 
field from the nonlinear displacement field can be nontriv- 
ial without further priors. In the middle panels, we see that 
MAK again agrees extremely well with SIM NL with a very 
small bias: the scatter around the diagonal is much smaller 
than in the right panels although we see, in the 3 ft -1 Mpc 
smoothing case, an increase of the dispersion for positive val- 
ues of — V q ■ v and a slight bias at larger values of —V, ■ v. 
This comparison between SIM NL and MAK confirms fully 
the nontrivial nonlinear nature of MAK reconstruction. The 
left panels are very similar to the right panels, as expected. 
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Figure 7. Same as Fig. [5] but for the projected negative divergence of the displacement field (proportional to the over-density in the 
linear regime). MAK reproduces very well the SIM NL field, at all densities. On the contrary, there are significant differences in regions 
of medium to high densities between the MAK and the SIM IC fields: v = 1.5 to 2 regions are more patchy in MAK and SIM NL than 
in SIM IC, and the sharp, highest density peaks in SIM IC are not reproduced in the other two maps. 
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Figure 8. Same as Fig. [3] but for the divergence of the displacement fields. The resolution is 3 (resp. 8) h~ 1 Mpc in the upper (resp. 
lower) panel. Again, MAK reproduces very well the SIM NL displacement field, even though the curves are more noisy than in Fig. [3 
the consequence of taking the derivatives. 
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Figure 9. Upper left panel: z-projected simulated spatial distribution of particles of the third most massive cluster (with total mass 
Mtot = 1.4 X lO 15 /! -1 Mq and 4316 particles) of the Gaussian simulation at the present time. In this Figure we use the full resolution 
analysis, i.e., sample So- Upper right panel: z-projected Lagrangian region in the uniform grid for the particles ending in the cluster and 
the Zel'dovich SIM IC velocities. Lower left panel shows the MAK-reconstructed Lagrangian region and MAK-reconstructed displacement 
field. Lower right panel shows the true Lagrangian region (as in the top right panel) but with the SIM NL displacement field. We use 
arbitrary units but a common scale for all plots, so that the amplitudes can be compared. In the right top and right bottom panels we 
use the exact initial positions provided by the simulation and in the left bottom panel the Lagrangian region reconstructed by MAK. 



4.3 Present-day peculiar velocity fields 

An important and direct outcome of MAK reconstruction, 
the peculiar velocity field vmak(x), is a surprisingly good 
approximation to the present day true peculiar velocity field 

Vpoc(x). 

vmak (x) for a given particle is obtained simply by scal- 
ing of the MAK displacement obtained previously at recon- 
structed particle's position (node of the Lagrangian grid), 
and assigning the displacement to this particle's position at 
z = following equations Q and JSJ. It is unclear, except 
maybe at very large scales, that vmak(x) can approximate 



Vpec(x). Fig. H2l shows that this is indeed the case on scales 
of 8 Mpc. There, we show the z = v x -v y velocity field 
extracted at z = from a 20 Mpc -thick slice normal to 
the z-direction and cutting through the centre of the Gaus- 
sian simulation. (We use the Si sample in this Figure.) The 
left (resp. right) panel shows the v pec (x) field given by the 
simulation (resp. vmak(x)), smoothed with a top-hat kernel 
of radius 8 h^ 1 Mpc. The underlying colour map gives the 
simulated dark matter density field smoothed on the same 
scale and then projected; it is the same for both panels. On 
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Figure 10. Same as Fig.|I]but for the divergence of the displacement field. The left and right figures correspond to 3 h^ 1 Mpc and 8 
h —1 Mpc smoothing scales respectively. Both MAK and SIM NL which agree very well have a cutoff at large value of v indicated by the 
vertical line (see text for further explanations). 
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Figure 11. Scatter plots comparing in a point- by-point perspective (nodes of a regular grid) the density of the MAK field to the density 
of the SIM IC field (left panels), the MAK vs. SIM NL density fields (middle) and the SIM NL vs. SIM IC (right). The upper (resp. 
lower) panels are for a 3 (resp. 12) h~ 1 Mpc smoothing length. We have used the negative divergence of the displacement as a proxy for 
the density. We recover the expected tight correlation between MAK and SIM NL and the absence of any significant bias. A horizontal 
line is displayed on each panel. It corresponds to the expected upper bound — V q • v = 3 as discussed in the text. In the bottom right 
panel, the curve shows the result given by the spherical collapse approximation [see Eq. 11 H i . The limitations of this model, discussed 
later in Section|H] justify the choice of 12h~ 1 Mpc instead of 8/i _1 Mpc for the smoothing length. 
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Density and velocities in simulation Density in simulation, velocities from MAK 




x [ Mpc ] x [ Mpc ] 

Gaussian z=0 

Figure 12.2 = v x -v y velocity field of a 20 h _1 Mpc -thick slice normal to the z-direction cutting at a depth of ~ 1/3 of the simulation 
box in z through the Gaussian simulation; the left and right panels respectively correspond to the simulated (peculiar) velocity field and 
to the reconstructed velocity field. The 8 h~ 1 Mpc smoothing scale guarantees that we are not affected by shell crossing. Underlying the 
velocity field we show in colour map the 3-dimensional density field at z = cut from the same z-slice and projected in z (it is the same 
density field in the left and right panels, and black is over-dense). There is quasi-perfect agreement between the two velocity fields. 



these scales, we find excellent agreement between the simu- 
lated peculiar velocities and the reconstructed velocities. 

On much smaller scales, Fig. 1131 details the simu- 
lated and reconstructed present-day (v x ,v y ) components of 
vmak(x) and v pec (x) for particles selected in a 4 3 ft -1 Mpc 
region of the Gaussian simulation. (We use the So sample 
here.) Although the agreement in the sparse/under-dense 
regions is still generally very good, differences arise between 
the simulated and the reconstructed particle velocities in 
regions of medium and higher densities. 

A thorough comparison between MAK and other meth- 
ods of proper-velocity evaluation (e.g. Yahil et al 1991 ; 
Valentine et al. 2000 and Shaya et al. 1995) remains to be 
done in future works where we plan to apply MAK to mock 
and large-scale redshift galaxy catalogues. 

The very good agreement on scales of ~ 8 ft -1 Mpc 
suggests a possible use of MAK to turn present-day redshift 
catalogues into real space catalogues, in an iterative process 
to deal with redshift space distortion. In applying MAK to 
real or mock galaxy catalogues, a quantitative error analysis 
is surely needed to access the goodness of MAK in evalua- 
tion of quantities such as cosmological parameters or mass 
to light ratio (some of these issues have been discussed in 
Mohayaee and Tully 2005 for the application of MAK to 
Near By Galaxies (NBG) catalog). 

4.4 Possible resolution issues 

Our study of the resolution effects consisted of analysing 
the full resolution sample (So) and dense sample (52) [see 
Section EH! and Table 1] and of comparing the results to the 



sparse sample (Si) used for most of the analyses up to now. 
We did not notice any significant differences between the full 
and sparse resolutions. 

To illustrate our point, the top and bottom panels of 
Fig. [Til compare the PDFs of SIM IC, MAK and SIM NL 
density fields smoothed with a top-hat kernel of radius 3 
and 8 ft -1 Mpc respectively using the full set So of 128 3 
reconstructed particles. There are no significant differences 
with respect to the PDF of Fig. HOI computed for the set Si- 
Note that this nice agreement between the sparse and full 
samplings is further supported by the right panel of Fig. 
and the last column of Table 1: the "ideally" reconstructed 
Lagrangian positions in Si are consistent with the fraction 
of the particles in So that are reconstructed with accuracy 
better than two mesh sizes. 

The dense sample, S2, tests as well, to some extent, the 
importance of tidal and edge effects on the reconstruction. 
In practice, this can have important consequences for recon- 
structions of volume-limited (in the literal sense) catalogues 
which might not extend well above the scale of homogene- 
ity. For the dense sample size considered here, lOO/i" 1 Mpc , 
these effects seem to be negligible. For instance, the fraction 
of "ideally" reconstructed particles is 17% to be compared 
with the 18% obtained for the full sample, So, reconstruc- 
tion. 

Note that we have performed these resolution tests as 
well for the x 2 model and the dense sample for the PVM 
and have obtained similar results. 
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x coordinate (Mpc/h) 



Figure 13. left panel: z = V^-Vy velocities of particles selected from a 4 h~ 1 Mpc -side and 20/i — 1 Mpc deep sub-volume taken from 
the Gaussian simulation. Simulated peculiar (resp. reconstructed) velocities are shown with the red (resp. black) arrows. There is fairly 
good agreement in the voids, but more noticeable differences in the outskirts of haloes, right panel: scatter plot of the reconstructed 
and simulation velocities at z=0 for the Gaussian sample of 64 3 , sample Gaussian Si of Table 1 , demonstrates closely the goodness of 
the reconstructed velocity field. Axes are in box units (one unit is 200 h _1 Mpc ). 
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Figure 14. Same as Fig. 1101 but for the full set So of 128 3 particles rather than its sparse sample, Si. 



5 RECONSTRUCTING NON-GAUSSIAN 
DENSITY FIELDS 

In this section, we examine different non-Gaussian initial 
conditions and verify if the results obtained previously for 
the Gaussian case still hold. We shall see that this is indeed 
the case, namely that the displacement fields reconstructed 
by MAK match extremely well their true nonlinear counter- 
parts given by the simulations. 

However, the real goal here is to study the properties 
of the initial density field or in other words the divergence 
of the initial displacement field which differs from the total 
final field, that is best reconstructed by MAK: the nonlinear 
contribution due to gravitational collapse should be sub- 
tracted from the reconstructed displacement field to allow 
tests for primordial non-Gaussianity. This will be confirmed 
by the subsequent analyses which will demonstrate that a 
"simple" use of MAK reconstruction is insufficient for re- 



covery of a small level of non-Gaussianity. By "simple" we 
mean the extrapolation of the nonlinear displacement field 
to early epochs using the Zel'dovich approximation to infer 
the initial density field. 

To carry out our analyses and reach these conclusions, 
we have explored the large spectrum of non-Gaussian mod- 
els, detailed in Section 



(i) The % 2 model, examined in Sec. 15. II which is strongly 
non-Gaussian but with the same topological properties as 
its Gaussian seed; 

(ii) The primordial voids model (PVM), studied in Sec. 
15.21 which is strongly non-Gaussian and initially inhomoge- 
neous; 

(iii) The weakly non-Gaussian Quadratic (Q) models, dis- 
cussed in Sec. 15.31 
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5.1 x 2 model 

In Fig. 1151 we examine —V ■ v in a 10 h^ 1 Mpc thick slice, 
similarly as in Fig. Q while Fig. Ilbl compares the PDF of 
-V-v for SIM IC, MAK and SIM NL. The measurements are 
performed using a top-hat smoothing length of 8 h~ l Mpc. 
The PDFs of MAK and SIM NL agree perfectly, which is 
consistent with the visual inspection of Fig. 1151 They differ 
only slightly from SIM IC PDF. This might be surprising at 
first sight, because visually SIM NL and MAK look rather 
different from SIM IC, at least for large values of —V ■ v. 
As explained in length in our study of the Gaussian model, 
this difference is mainly due to the presence of collapsed 
objects which induce an abrupt cut-off in the tail of the 
PDF for large values of —V ■ v (This cut-off lies outside the 
plot range in Fig, ll6[l . However, a more careful inspection of 
Fig. 1151 shows a much better agreement in regions of lower 
density, which explains the very good match between all of 
the three PDFs. 

Except for the abrupt cut-off, the small difference as 
compared to the Gaussian case between the MAK/SIM NL 
PDF and SIM IC is due to two factors. First, the contribu- 
tion from initial non-Gaussianity is strong enough to com- 
pletely dominate over the contribution from gravity. Second, 
in this model the effect of gravitational clustering is less in- 
tense than in the Gaussian model because of the lower value 
of erg. A noticeable consequence of this lower normalisation 
is that there are less shell crossings and therefore the frac- 
tion of "ideally" reconstructed particles is much higher than 
for the Gaussian model (see last column of Table 1). 

5.2 PVM model 

In Fig. 1171 we examine —V ■ v, in a 10 Mpc thick slice, 
similarly as in Fig. Q while Fig. 1181 compares the PDF of 
-V • v for SIM IC, MAK and SIM NL. The measurements 
are performed with a top-hat smoothing length of 8 ft -1 Mpc 

Recall that for PVM, initial conditions are strongly in- 
homogeneous at very high redshifts. As a result, it is in 
principle meaningless to use MAK reconstruction. However, 
at some point, the initial conditions had to be almost homo- 
geneous, prior to the seeding of the primordial voids. Hence, 
we can still assume that we are starting from Lagrangian po- 
sitions setup on a homogeneous grid. In our case, this homo- 
geneous grid is the one where particles are located prior to 
the Zel'dovich displacement and void creation as explained 
in Section \Z. 31 With this in mind, a fair comparison between 
the middle and right panel of Fig. ll7l can be made and again 
the visual agreement between MAK and SIM NL is excellent. 
Expectedly, due to the early strong nonlinearities induced 
by voids creation and subsequent shell crossings, the frac- 
tion of "ideally" reconstructed particles is rather low (28%, 
see Table 1). 

From the previous arguments, comparing MAK/SIM 
NL to SIM IC has to be done with care. Clearly the strong 
"density contrasts" observed around voids in the left panel 
are expected to be smeared out in the middle and right pan- 
els. Keeping this in mind, we notice that there is a global 
large scale qualitative agreement between the three panels 
of Fig. 1171 including the positions and extensions of the 
under-dense regions, although the amplitudes of the fluctu- 



ations are different. This means that our nonlinear displace- 
ments fields have kept at least some informations on the 
non- Gaussian nature of initial conditions even if the middle 
and right panels differ significantly from the left one. 

The results of the previous discussion are nicely illus- 
trated by Fig. EHI that we can now comment on in detail. 
The bell-shaped part of the SIM IC PDF corresponds to 
the initial Gaussian field generated prior to the seedings of 
the voids. The strong over-dense and under-dense regions 
appearing during the creation of the voids increase the vari- 
ance of the density field by adding tails to the PDF both in 
the high and the low parts. This explains why the bell-shape 
of the SIM IC PDF on Fig. ll8l is in fact narrower than a pure 
Gaussian with the same variance, given the choice of repre- 
sentation, in units of — V 9 • v/cr. Note that this additional 
contribution to the variance is, at least in the linear regime, 
a pure transient, as discussed furthermore at the end of Sec- 
tion |S] If one neglects additional non-Gaussianity brought 
about by gravitational instability, one thus expects that the 
bell-shaped part of the PDF converges to the pure Gaussian 
at late times, as one can indeed observe in Fig. 1181 

The MAK and SIM NL PDFs agree with each other 
very well as expected. In the positive — V 9 • v part, nonlin- 
ear gravitational clustering and subsequent shell crossings 
destroy information related to the strong over-dense spheres 
around the voids. As a result, in this regime the shape of the 
PDF is qualitatively very similar to what would be obtained 
from the Gaussian case. However, as noticed while examin- 
ing Fig. 1171 the informations on the primordial voids is still 
present in the tail for the low values of — V g ■ v. Clearly, de- 
spite the fact that SIM NL/MAK PDFs differ enormously 
from SIM IC PDF, they have still preserved a detectable 
signature of primordial non-Gaussianity. 



5.3 Quadratic models 

Fig.EUshows the PDF of -V ■ v for the SIM IC, MAK and 
SIM NL, with a top-hat smoothing of 8 h~ x Mpc for the 
series of quadratic models Q (recall that these are models 
with very small non-Gaussianity). These models evolve to 
very similar nonlinear stages since they are all weakly non- 
Gaussian: gravity dominates significantly over initial non- 
Gaussianity suggesting that it will be nontrivial to recover 
the latter. Once again, the reconstructed nonlinear displace- 
ment field matches the simulated one for all of these models. 
As expected, the fraction of "ideally" reconstructed particles 
is similar to that in the Gaussian model (41 — 43% versus 
37%, see last column of Table 1). 



6 SEPARATION OF THE INITIAL AND 
GRAVITATIONAL-INDUCED 
NON-GAUSSIANITIES: THE SPHERICAL 
COLLAPSE MODEL 

We see that MAK reconstructs the nonlinear displacement 
field with a tremendous accuracy. The problem is that, in 
order to constrain the primordial non-Gaussianity, it is es- 
sential to recover the initial density contrast, i.e. the diver- 
gence of the linear displacement field. To achieve this goal 
it is necessary to model the effect of gravitational instabil- 
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Figure 15. Same as lower panel of Fig. [7] but for the x 2 model. 
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Figure 16. Same as Fig. HUI but for the \ 2 model. 



ity on the displacement field. Several approaches might be 
chosen: 

(i) Semi- analytic field reconstruction method in an it- 
erative way. The method would inspire closely from the 
ideas developped in the building of ZTRACE reconstruc- 
tion (Monaco & Efst athiou 1999), i.e. wo uld be compa- 
rable to Peak-Patch jBond fc Mversl 1 199611 or Pinocchio 
jMonaco. Theuns fc Taffonil 12002 ) , but in a reverse way. 
More explicitely, the idea is to use a realistic approximation 
of the dynamics depending in a unique way on the initial 
displacement field and the goal is to reconstruct the latter 
in an iterative manner. For instance, one can use Lagrangian 
perturbation theory up to some order, the latter being de- 
termined by the level of non-Gaussianity one aims to probe, 
e.g. the skewness (second order needed) or the skewness and 
the kurtosis (third order needed) of the initial distribution 
function. The subtlety of this approach is that one needs to 
truncate the sought-after initial fluctuations at some scale in 



order to deal appropriately with shell crossings in over-dense 
regions. The corresponding smoothing scale depends on lo- 
cation, according to the mass scale of the collapsed object 
formed at present time. 

If one uses Lagrangian perturbation theory, the nonlinear 
displacement field is indeed, for the growing mode, entirely 
determined as a function of the initial one, even if one has 
to take into account the adaptive smoothing of initial condi- 
tions needed for dealing with collapsed regions, as has just 
been explained. In this approach, the goal would then be to 
solve an implicit equation on the initial displacement field by 
equating the subsequent theoretical nonlinear displacement 
field with the reconstructed one by MAK. To solve this im- 
plicit equation, the most obvious approach would rely on an 
iterative method, which is unfortunately not guaranteed to 
converge. 

(ii) Statistical method relying on loop perturbation theory. 
A less elegant, but simpler and still powerful method con- 
sists of working directly on the PDF of the divergence of 



© 0000 RAS, MNRAS 000, 000-000 



Reconstruction of primordial density fields 19 





Figure 18. Same as Fig. I16l but for the PVM model. 



the reconstructed displacement field to infer the PDF of 
the initial one. Again, some modelling of the dynamics is 
needed, to infer the mapping between the former and the 
latter. The best way to achieve that relies on perturba- 
tion theory and its one-loop (or h igher order) corrections 
fe.g.. IScoccimarro fc Friemanlll996l) . The advantage of this 
approach is that it allows one to probe scales which are close 
to the nonlinear regime. Since MAK reconstruction is very 
good, even at mildly nonlinear scales, this method would be 
optimal, but still quite involved from the analytic point of 
view. 

(iii) Statistical method relying on the top-hat spherical 
model. To simplify furthermore the approach, with still 
in mind the goal of remapping of the probability dis- 
tribution, one can use the top -hat spherical model (e.g., 
iFosalba fc Gaztanaeal Il998al lbf) . This approximation has 



been demonstrated to work very well, both from the theoret- 


ical (Bernardeau 1992, 1994a 




iGaztanaea & Fosalball998l 


IFosalba & Gaztanaeal Il998allbl 


and the numerical points 



of view iFosalba fc Gaztanaealll998allbl . However, it works 
only in the regime where fluctuations of the considered field 
(here the divergence of the displacement) are very small. In- 
deed, this model presents wrong loop corrections, but has 
the considerable advantage of relying on a very simple for- 
malism. 

It would be beyond the scope of this paper to try to apply 
approaches (i) and (ii) to our data, because they are both 
rather involved, from the numerical and analytical points 
of view. Instead, we concentrate on the spherical top-hat 
model, to demonstrate that it is possible to estimate the 
level of non-Gaussianity of the initial displacement field from 
the MAK-reconstructed one. 

Under the top-hat spherical approximation, the diver- 
gence of the displacement field can be expressed as a function 
of the density contrast as follows: 

i> = V,v = 3[(l + 5)- 1/3 -l] . (9) 
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Figure 19. Same as Fig. 1101 but for the series of quadratic models. Q-ioOi Q+iOOi Q— 50i Q+50 respectively correspond to the upper 
left, upper right, lower left and lower right panels. 



We now relate tp to the initial linear density contrast, ipL = 
—Sl- In the spherical top-hat framework, the approximation 



1 + 8: 



-3/2 



(10) 



relates 8 to S L iBernardeaulll994bl) . In the above equation, 
the time dependence of the growing mode extrapolated to 
present time has been included in Sl- This approximation 
turns out to be excellent, independently of t he value of the 
cosmological parameters feernardeaulll994bl) . 
We thus obtain 

1/2 

- 1 



1> 



1 + 



(11) 



In the spherical approximation framework, this gives us, in 
Lagrangian space, the transformation from the initial diver- 
gence of the displacement field to the final one. This formula 
is unchanged if top-hat smoothing is appli ed to the fie lds, 
here tpL and tp, under consideration (IBernardeaul l 1994a ) . 
Inverting Eq. 11 U . we obtain the simple relation 



/ 3 
^= 2 



1+ ! 



i 



(12) 



The spherical collapse approximation has already been 
used to explicitly calculate the PDF in a local Lagrangian 



formalism (Protogeros & Scherrer 1997; Protogeros et al. 
1997). The above equation gives us the change of variable to 
obtain the PDF of the initial divergence of the displacement 
field as a function of the non-linear, reconstructed one, by 
simply using 



P(ip L )d^ L = P(tp)dtp. 



(13) 



Note that the relation between tp and ipL [Eq. 11H 1 imposes 
an upper bound on the possible values of tpL, tpL ^ —3/2, 
which in turns implies that tp ^ —3, as already discussed 
extensively in Section 14.21 This also shows the limits of 
the spherical collapse approximation: it is expected to fail 
for collapsed objects, i.e. when tp approaches the bound 
tp = —3, as illustrated by Fig. 1111 There is in fact cer- 
tainly a way to fix this by using other approximations in 
the highly nonlinear regime that we do not explore here 
(e .g. using the so called "h alo model" for instance discussed 
in IScoccimarro et al.ll200lT) . 

Here we test the spherical collapse model approxima- 
tion applied to the divergence of the MAK-reconstructed 
displacement field at a scale of 12/i _1 Mpc, where we know 
that one-loop corrections are expected to be negligible 
and thus where the spherical collapse model should per- 
form well, even in the non- Gaussian case, as argued by 
iFosalba fc Gaztanaeal Jl998rJ) . 
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We first examine in Fig. 1201 the Gaussian, x 2 and Q 
models. We shall pay attention to the special case of the 
PVM model at the end of this section. The success of the 
spherical model is unquestionable. Except in the tails (es- 
pecially the right-handed one, as argued above), the shape 
of the initial PDF is recovered to very good accuracy. As 
discussed extensively in the caption of Fig. 1201 even subtle 
features such as those observed in the Q models are repro- 
duced, lending credence to the fact that it seems possible 
to detect even small levels of primordial non-Gaussianity, 
at least with MAK reconstruction. Of course, this assertion 
does not take into account the instrumental and observa- 
tional effects as shall be discussed later in the conclusion. 
Furthermore, our approach here is rather phenomenological, 
since we do not present a serious treatment of the statistical 
significance of the comparison. 

We now turn to the PVM model. This model deserves 
a special treatment because the initial velocity field is not 
set up according to the Zel'dovich approximation, as ex- 
plained in detail in Section f2.3l As a consequence, there are 
transients even for the variance of the displacement field: 
in the linear regime, the expanding voids are in fact pure 
transients, but they induce strong density contrasts that in- 
crease the variance of the fluctuations, thus the variance of 
the divergence of the displacement field. A consequence is 
that standard linear theory, with only the growing mode, 
does not apply. For instance, the value of a (where a 2 is the 
variance of the field) at 12/i _1 Mpc measured in the initial 
—V ■ v scaled to present time is 0.94, larger than the mea- 
sured one from the nonlinear displacement field, 0.68. The 
spherical collapse model captures only the growing modes: 
they are determined by the initial Gaussian field prior to the 
void seeding. This one had a value of a of 0.71, hence of the 
same order as the present one, as expected. In Fig. 1211 we 
thus assumed for the comparison between true initial — V ■ v 
PDF and predicted one a value of a of 0.71 for the latter. 
In practice, note that the test for non-Gaussianity does not 
need the precise knowledge of this initial a: only the shape 
of the PDF matters, and has to be compared with the best 
fitting Gaussian in a trustworthy range. 

Here the match between the predicted and the initial 
PDF is not so good anymore, although some non-Gaussian 
features seem to be captured, namely the right hand ex- 
cess compared to the Gaussian case, but this might as 
well be a mere coincidence. The problem here is that the 
arguments in favour of the sp herical top-hat model (e.g., 
iFosalba fc Gaztanagalll998allrl might not apply to the PVM 
model and its very strong transient nature, despite the fact 
that the voids are chosen here to be perfectly spherical. 2 
The validity of Eq. [11 1 1 and especially of the assertion that 
smoothing does not affect this equation is certainly ques- 
tionable for the PVM model. 



7 DISCUSSION 

In this work, we have tested the Monge-Ampere-Kantorovich 
(MAK) reconstruction against JV-body simulations. We have 

2 but they fragment in a nontrivial way due to the additional 
fluctuations induced by the Gaussian field. 



examined the standard ACDM model and 6 additional mod- 
els with non-Gaussian initial conditions: a \ 2 model where 
the initial density field is the square of a Gaussian, a model 
with primordial voids, and four mildly non-Gaussian mod- 
els where the primordial gravitational potentials includes a 
small, quadratic correction. 

The main results of this paper are as follows: 

• In its essence, MAK is supposed to reconstruct the non- 
linear displacement field between initial and final positions. 
The analyses of this paper show that it achieves this with an 
unprecedentedly high degree of accuracy, at least at scales 
larger than the size of rich clusters, i.e., scales above a few 
Mpc. In particular, it captures in a nontrivial way the non- 
linear contribution from gravitational instability, well be- 
yond the Zel'dovich approximation. 

• From the MAK-reconstructed nonlinear displacement 
field, one can infer the present velocity field, using the 
Zel'dovich approximation. Again, our iV-body simulation 
analyses demonstrate the success of MAK in fulfilling this 
goal. 

• The displacement field provided by MAK can be used 
as well to constrain the statistical nature of the primor- 
dial density field. But, since MAK encodes nonlinearities 
due to gravitational clustering, it is difficult to disentan- 
gle these dynamical contributions from possible initial non- 
Gaussianities. However, as illustrated here in a simple way 
with the spherical collapse model, we show that it is possible 
to truly recover the statistical properties of the primordial, 
linear density field using additional modeling of the dynam- 
ics. Furthermore, we envisage possible improvements to be 
made by elaboration of MAK to a more sophisticated mod- 
eling of the dynamics combined with adaptive smoothing of 
the initial fluctuations [similarly as in ZTRACE (Monaco & 
Efstathiou 1999)]. 

The results presented in this paper are obtained in 
the "idealized" framework of iV-body simulations. Extra 
complications arise in reconstruction from real galaxy cat- 
alogu es (for application of MAK t o real galaxy catalogues 
see IMohavaee fc Tullvl 120041. EobH) . Here, we discuss red- 
shift space distortion, edge-effects, biasing and catalog in- 
completeness, which are the most relevant issues for MAK 
reconstruction: 

(i) Redshift space distortion: MAK has already been di- 
rectly applied to redshift catalogues by modifying the cost 
function © using once again the Zel'dovic h approxima- 
tion (for detail, see IMohavaee fc Tullvl 120041) . An alterna- 
tive method would be to deduce real space positions from 
redshift space data, using MAK in an iterative way. Indeed, 
we have shown in paragraph 14.31 very good agreements be- 
tween the present-day simulated peculiar velocities and the 
MAK-reconstructed ones already on ~ 8 h -1 Mpc scales. 

(ii) Boundary effects: in this paper, we have partly ad- 
dressed edge and tidal effects using the dense samples. We 
have noticed that these effects are small but the volume con- 
sidered was large, a cube 100 h~ 1 Mpc aside. Clearly, large 
scale tidal effects and wrong assignments near the edges of a 
catalog would be important if it had a small volume coverage 
or/and intricate boundaries. 

(iii) Biasing: in a real galaxy catalog, there is the prob- 
lem of biasing that is the relationship between the present 
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Figure 20. Comparison between the PDF of the true initial — V- v (SIM IC, dots) to the one inferred from the reconstructed displacement 
field using the top-hat spherical collapse approximation (SIM IC, predicted, thick solid curve), after a smoothing with a top-hat window 
of radius 12h~ 1 Mpc. The choice of 12h~ 1 Mpc instead of 8h~ 1 Mpc is justified by the plots of Fig. 1111 There are additional curves as 
detailed in the panels: the PDFs of — Vq ■ v for the nonlinear displacement fields, measured in the simulation (SIM NL, dot dashes) and 
reconstructed (MAK, dashes), and finally the pure Gaussian prediction (Gaussian, analytic, thin solid curve, not present on the upper 
right panel). Each panel corresponds to a given model as indicated on the figures. 

The important curves to examine on each panel are the dotted ones (initial PDF) and the thick solid ones (predicted initial PDF from the 
MAK nonlinear PDF), while keeping the Gaussian limit in mind: they should in principle superpose if the combination MAK+spherical 
collapse model works. This is indeed the case in all panels, except in the tails, especially the right-hand side one, but only for rather 
large values of | — V 9 • v/iry-v] ~ 2. Indeed, in thi s regime, especially for positive values of — V 9 • v, the spherical collapse model is 
known to have limitations (e.g.. lBernardeaulll994bTl . Except for this, the agreement between the prediction and the measurement is 
excellent. It is in fact so good, that for the Q models, the weakly non-Gaussian features are recovered, if one examines the upper part of 
the curves displayed in the four lower panels. The curves displayed on each are however slightly irregular, because the scale considered 
here is a significant fraction of the box size and the measured PDFs are thus expected to be contaminated by finite volume effects. This 
explains for instance the fact that in the Gaussian case, the dotted curve does not exactly superposed to the thin solid curve, i.e. the 
pure Gaussian prediction. Note that even this small deviation is reproduced by the MAK+top-hat prediction, at least for the right-hand 
part of the bell-shaped curves. 
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total matter distribution and the present observed light dis- 
tribution. For each observed galaxy of a given luminosity, L, 
one has to infer a given mass which for instance can be a 
simple function of L. In this way, one can transform a dis- 
tribution of points with given luminosities to a distribution 
of points with given masses. A supplementary complication 
arises since MAK requires in essence all the points to have 
the same mass. To achieve this, one can break the galax- 
ies into equal mass particles, place them at the position of 
the galaxy and perform MAK reconstruction. 3 In the hier- 
archical clustering framework, galaxies are located in dark 
matter halos, that is in condensed objects. We know from 
the results of this paper that MAK reconstructs very well 
the nonlinear displacement field even for regions that have 
experienced shell-crossings. In the latter case, the only prob- 
lem is the wrong assignment of the particles inside the col- 
lapsed objects, but this does not affect the reconstruction 
of the displacement field significantly, and hence the veloc- 
ity field inferred from it using the Zel'dovich approximation. 
Therefore, provided that the right mass has been given to 
each galaxy and that the catalog is complete enough so that 
the underlying total mass distribution is well-traced, the re- 
constructio n is expected to perform well iMohavaee fc Tullvl 
|2004 120051) . There is however a subtlety which has to be 
taken into account, during mass assignment to galaxies: a 
different treatment has to be performed for field galaxies 
which have their own halo and galaxies belonging to rich 
clusters. In the latter case, what matters is the mass of the 
halo of a cluste r and not the mass of the "sub-halos" of each 
of its galaxies dMohavaee fc TullvlEoOsT) . 

(iv) Luminosity segregation: the last issue to consider 
is catalog incompleteness. For instance, in a standard 
magnitude-limited catalog, there is a luminosity segregation 
effect which arises due to fewer and brighter objects with in- 
creasing distance from the observer. It is possible to compen- 
sate ^fo^t_bv_j ; s^i^niiu^ al_arger mass to more distant galax- 
ies ijMohavaee fc Tullvll2005Fl : this means that two galaxies 
of the same luminosity but at different distances from the 



3 In practice, one distributes the particles into Gaussian clouds 
with negligible dispersion around the original galaxy. One reason 
for this is to speed up the assignment algorithm. 



observer will be given different masses. With this procedure, 
the reconstruction is in principle expected to be as least bi- 
ased as possible. However, clearly the signal-to-noise ratio 
would decrease with distance from the observer, and more 
subtlely, the expected strong dependence of galaxy cluster- 
ing with luminosity might complicate the interpretation of 
the results: obviously, it is needed to have catalogs complete 
enough in the faint end of the luminosity function in order to 
treat appropriately the biasing issues discussed in previous 
point. 
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